# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 08:18:36 2022

@author: Derek Rodriguez Pacheco, PhD 
IUSS Pavia
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
import matplotlib as mpl
from Spectra.SpectrumA import spectrumAD 

mpl.rc('font',family = 'Times New Roman')

m= 950                           ## Mass in N*s2/m or kg
g= 9.81                          ## Acceleration of gravity in m/s2

adjfactors= np.array(pd.read_csv('Adjustment_Factors_Displacement_CM.txt', header= None, delim_whitespace= True)[0])
indexesb= np.array(pd.read_csv('Indexes_Beginning.txt', header= None, delim_whitespace= True)[0])
indexesf= np.array(pd.read_csv('Indexes_End.txt', header= None, delim_whitespace= True)[0])

## Model's parameters

k1= 0.5
k2= 0.01
sigAct= 3.8
betaf= 0.51
epsSlip= 0
epsBear= 0
rBear= 1


# s1p= 3*afs
s1p= 2
# e1p= 2*afd
e1p= 6
# s2p= 9*afs
s2p= 4
# e2p= 10*afd
e2p= 25
# s3p= 35*afs
s3p= 15
# e3p= 55*afd
e3p= 50
pinchX= 0.5
pinchY= 0.01
damage1= 0
damage2= 0
beta= 0

s1n= s1p*31
e1n= np.copy(e1p)
s2n= s2p*31
e2n= np.copy(e2p)
s3n= s3p*31
e3n= np.copy(e3p)

width= 600
height= 1053

testnumber= [13,14,15,16,17,18,19,20,24,25,26,27,28,40,41]

##

datacm= {}
# acceltcm01= {}
# acceltcm02= {}
accelcm= {}
timecm= {}
timecmc= {}

data= {}
datad= {}
accelda= {}
dataR = {}
# accelda1= {}
# accelro= {}
accelbase= {}
dispda= {}
dispdacm= {}
timearray= {}
energyarray= {}
energyarraycm= {}
maxaccel= []
minaccel= []
maxaccelf= []
matchindexmin= []
matchindexmax= []
matchdispmin= []
matchdispmax= []
timetotal= []
acceltotal= []
dispnet= {}
dispnetcm= {}
disp1 = {}
disp2 = {}
baserot = {}
# timestart= []

for i in range(len(testnumber)):
    datacm['%d' % (i+1)]= np.array(pd.read_csv('Data_Center_Mass\\accelAndDispl_15Hz\\%d_AccelAndDispl_15Hz.txt' % testnumber[i], header= None, delim_whitespace= True))
    accelcm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,4])[indexesb[i]:indexesf[i]+1]
    timecm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,0])[indexesb[i]:indexesf[i]+1]
    timecmc['%d' % (i+1)]= timecm['%d' % (i+1)]-timecm['%d' % (i+1)][0]
    disp1['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,5])[indexesb[i]:indexesf[i]+1]
    disp2['%d' % (i+1)] = np.array(datacm['%d' % (i+1)][:,6])[indexesb[i]:indexesf[i]+1]
    
    
for i in range(len(testnumber)):
    dm = []
 
    for j in range(len(disp1['%d' % (i+1)])):
        d1 = disp1['%d' % (i+1)][j]
        d2 = disp2['%d' % (i+1)][j]
        dm.append((d1+d2)/2)
        
    dispdacm['%d' % (i+1)] = np.array(dm)*adjfactors[i]


data['1']= np.array(pd.read_csv('13_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['2']= np.array(pd.read_csv('14_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['3']= np.array(pd.read_csv('15_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['4']= np.array(pd.read_csv('16_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['5']= np.array(pd.read_csv('17_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['6']= np.array(pd.read_csv('18_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['7']= np.array(pd.read_csv('19_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['8']= np.array(pd.read_csv('20_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['9']= np.array(pd.read_csv('24_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['10']= np.array(pd.read_csv('25_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['11']= np.array(pd.read_csv('26_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['12']= np.array(pd.read_csv('27_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['13']= np.array(pd.read_csv('28_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['14']= np.array(pd.read_csv('40_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['15']= np.array(pd.read_csv('41_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))

dataR['1']= np.array(pd.read_csv('13_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['2']= np.array(pd.read_csv('14_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['3']= np.array(pd.read_csv('15_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['4']= np.array(pd.read_csv('16_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['5']= np.array(pd.read_csv('17_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['6']= np.array(pd.read_csv('18_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['7']= np.array(pd.read_csv('19_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['8']= np.array(pd.read_csv('20_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['9']= np.array(pd.read_csv('24_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['10']= np.array(pd.read_csv('25_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['11']= np.array(pd.read_csv('26_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['12']= np.array(pd.read_csv('27_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['13']= np.array(pd.read_csv('28_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['14']= np.array(pd.read_csv('40_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
dataR['15']= np.array(pd.read_csv('41_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))


for i in range(15):
    accelbase['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,1])
    timearray['%d' % (i+1)]= data['%d' % (i+1)][:,0]-data['%d' % (i+1)][0,0]
    energyarraycm['%d' % (i+1)]= [np.trapz(accelcm['%d' % (i+1)][:j]*m*(-1),dispdacm['%d' % (i+1)][:j]/1000) for j in range(len(accelcm['%d' % (i+1)]))]
    baserot['%d' % (i+1)] = dataR['%d' % (i+1)][:,1]
    
cc= 0

for i in range(len(testnumber)):
    for j in range(len(timearray['%d' % (i+1)])):
        if i==0:
            timetotal.append(timearray['%d' % (i+1)][j])
        else:
            timetotal.append(cc+timearray['%d' % (i+1)][j])
        acceltotal.append(accelbase['%d' % (i+1)][j])
    cc= timetotal[-1]
    if i!=len(testnumber)-1:
        for j in range(int(round(15/timearray['%d' % (i+1)][1]))):
            if i==0:
                timetotal.append(timearray['%d' % (i+1)][-1]+(j+1)*timearray['%d' % (i+1)][1])
            else:
                timetotal.append(cc+(j+1)*timearray['%d' % (i+1)][1])
            acceltotal.append(0)
        cc= timetotal[-1]

#### SDOF Analysis

# for i in range(15):
np.savetxt('TH/Time_Total.txt', timetotal)
np.savetxt('TH/Acceleration_Total.txt', acceltotal)

recordsa= 'TH/Acceleration_Total.txt'
recordst= 'TH/Time_Total.txt'

dts= 0.001953
mass= m/1e6

dur= timetotal[-1]

fh= open('FMdata.tcl', 'w')

fh.write(
    'set FMfilet "%s";\nset FMfilea "%s";\nset dtg %f;\nset durs %f;'
    '\nset M %f;'
    '\nset s1p %f;\nset e1p %f;\nset s2p %f;\nset e2p %f;\nset s3p %f;'
    '\nset e3p %f;\nset pinchX %f;\nset pinchY %f;\nset damage1 %f;'
    '\nset damage2 %f;\nset beta %f;' 
    '\nset s1n %f;\nset e1n %f;\nset s2n %f;\nset e2n %f;\nset s3n %f;\nset e3n %f;'
    '\nset width %d;\nset height %d;'
    '\nset k1 %f;\nset k2 %f;\nset sigAct %f;\nset betaf %f;\nset epsSlip %f;\nset epsBear %f;\nset rBear %f;' % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height, k1, k2, sigAct, betaf, epsSlip, epsBear, rBear))
    # % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height))
fh.close()
subprocess.call('OpenSees MDOF_CalibrationMH_Grav.tcl', shell=True)

####

sdofd= {}
sdofa= {}
sdoft= {}
sdoff= {}
sdofe= {}
sdofr= {}

dtt = 0.001953

indexes_Beg = [0,int(42/dtt),int(85/dtt),int(127/dtt),int(170/dtt),int(212/dtt),int(255/dtt),int(298/dtt),int(342/dtt),int(387/dtt),int(431/dtt),int(475/dtt),int(518/dtt),int(562/dtt),int(606/dtt)]
indexes_End = [int(28/dtt),int(71/dtt),int(113/dtt),int(156/dtt),int(198/dtt),int(241/dtt),int(284/dtt),int(328/dtt),int(372/dtt),int(416/dtt),int(460/dtt),int(503/dtt),int(547/dtt),int(591/dtt),324987]
timei= [0,42,85,127,170,212,255,298,342,387,431,475,518,562,606]

for i in range(len(testnumber)):
    sdofd['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/DispT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])
    sdofa['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])
    sdoft['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[0][indexes_Beg[i]:indexes_End[i]])-timei[i]
    sdoff['%d' % (i+1)]= sdofa['%d' % (i+1)]*0.001*m*-1
    sdofe['%d' % (i+1)]= [np.trapz(sdoff['%d' % (i+1)][:j],sdofd['%d' % (i+1)][:j]*0.001) for j in range(len(sdoff['%d' % (i+1)]))]
    upR = np.array(pd.read_csv('ResultsAtt/UpliftT.out', header= None, delim_whitespace= True)[2][indexes_Beg[i]:indexes_End[i]])
    upL = np.array(pd.read_csv('ResultsAtt/UpliftT.out', header= None, delim_whitespace= True)[1][indexes_Beg[i]:indexes_End[i]])
    sdofr['%d' % (i+1)] = ((upL)-(upR))/width

'''    
lN = len(sdofd['%d' % 14])
lE = len(dispdacm['%d' % 14])

if(lN>lE):
    l = lE
else:
    l = lN
    
relerr = abs(np.array(sdofd['%d' % 14][3:l])-np.array(dispdacm['%d' % 14][3:l]))/abs(np.array(dispdacm['%d' % 14][3:l]))
re = np.mean(relerr[3:])   
'''
## Plotting

testnum2= [2,3,4,5,6,7,8,9,11,12,13,14,15,17,18]
intensity= [0.17,0.17,0.19,0.20,0.20,0.19,0.20,0.40,0.29,0.29,0.34,0.39,0.46,0.51,0.52]
location1= [2500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500]
location2= [2800,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500,7500]
location3= [-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58]
location4= [3.7,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2,5.2]
subplotnum = ['a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p']

nr= 3
nc= 5

fig1, ax1 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        ax1[i,j].plot(timecmc['%d' % count], energyarraycm['%d' % count], color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax1[i,j].plot(sdoft['%d' % count], sdofe['%d' % count], color= '#FF8787', lw= 1.5, label= 'Numerical model')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax1[i,j].set_xlabel('Time (s)')
            # ax1[i,j].set_xticks([0,10,20,30])
            # ax1[i,j].set_xlim([0,30])
        if count==1 or count==6 or count==11:
            ax1[i,j].set_ylabel('Absorbed Energy (J)')
            # ax1[i,j].set_yticks([0,800,1600,2400,3200])
            # ax1[i,j].set_ylim([0,3200])
        ax1[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, location1[count-1]), 
                      xycoords = 'data',
                      xytext = (+0, -5), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax1[i,j].set_xticks([0,10,20,30,40])
        ax1[i,j].set_xlim([0,40])
        
        ax1[i,j].annotate('%s)' % subplotnum[count-1], xy=(35,3700),fontsize=14)
        
        if count==1 or count==2 or count==3 or count==4 or count==5 or count==6 or count==7:
            ax1[i,j].set_yticks([0,1000,2000,3000,4000])
            ax1[i,j].set_ylim([0,4000])
            #ax1[i,j].set_yticks([0,100,200,300,400,500,600])
            #ax1[i,j].set_ylim([0,600])
        else:
            ax1[i,j].set_yticks([0,1000,2000,3000,4000])
            ax1[i,j].set_ylim([0,4000])
        
        ax1[i,j].set_facecolor('#DEE1E6')
        ax1[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax1[i,j].legend(loc=2)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig1.tight_layout()        
#fig1.savefig('MDOF_Absorbed_Energy.png', bbox_inches= 'tight',dpi=300)
fig1.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure10.tiff', bbox_inches= 'tight',dpi=300)





fig2, ax2 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        
        paccE = np.argmax(abs(accelcm['%d' % count]))
        paccN = np.argmax(abs(sdoff['%d' % count]))
        
        ax2[i,j].plot(dispdacm['%d' % count]*(-1), accelcm['%d' % count]*m, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax2[i,j].plot(sdofd['%d' % count], sdoff['%d' % count], color= '#FF8787', lw= 1.5, label= 'Numerical model')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax2[i,j].set_xlabel('Displacement (mm)')
        if count==1 or count==6 or count==11:
            ax2[i,j].set_ylabel('Force (N)')
        ax2[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (-95, location2[count-1]), 
                      xycoords = 'data',
                      xytext = (+0, -5), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        if(count==1):
            y = -86
        else:
            y = -120
        ax2[i,j].annotate('Relative error of peak = '+str(round(((abs(sdoff['%d' % count][paccN])-abs(m*accelcm['%d' % count][paccE]))/abs(m*accelcm['%d' % count][paccE]))*100,1))+'%',
                      xy = (-95, location2[count-1]-1000), 
                      xycoords = 'data',
                      xytext = (+0, y), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax2[i,j].set_xticks([-100,-75,-50,-25,0,25,50,75,100])
        ax2[i,j].set_yticks([-10000,-7500,-5000,-2500,0,2500,5000,7500,10000])
        ax2[i,j].set_xlim([-100,100])
        ax2[i,j].set_ylim([-10000,10000])
        
        ax2[i,j].annotate('%s)' % subplotnum[count-1], xy=(75,8000),fontsize=14)
        
        ax2[i,j].set_facecolor('#DEE1E6')
        ax2[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax2[i,j].legend(loc=2)
        
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig2.tight_layout()
#fig2.savefig('MDOF_Hysteretic_Response.png', bbox_inches= 'tight',dpi=300)        
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig2.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure12.tiff', bbox_inches= 'tight',dpi=300)

fig3, ax3 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        pcrdE = np.argmax(abs(dispdacm['%d' % count]))
        pcrdN = np.argmax(abs(sdofd['%d' % count]))
        
        ax3[i,j].plot(timecmc['%d' % count], dispdacm['%d' % count], alpha= 1, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax3[i,j].plot(sdoft['%d' % count], sdofd['%d' % count], alpha= 0.75, color= '#FF8787', lw= 1.5, label= 'Numerical model')
        #ax3[i,j].plot(timecmc['%d' % count][pcrdE],dispdacm['%d' % count][pcrdE],'o',mec='limegreen',mfc='none')
        #ax3[i,j].plot(sdoft['%d' % count][pcrdN],sdofd['%d' % count][pcrdN],'o',mec='#FF8787',mfc='none')
        
        
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax3[i,j].set_xlabel('Time (s)')
        if count==1 or count==6 or count==11:
            ax3[i,j].set_ylabel('Relative Displacement (mm)')
        ax3[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, -95), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax3[i,j].annotate('Relative error of peak = '+str(round(((abs(sdofd['%d' % count][pcrdN])-abs(dispdacm['%d' % count][pcrdE]))/abs(dispdacm['%d' % count][pcrdE]))*100,1))+'%',
                      xy = (1, -110), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax3[i,j].set_xticks([0,10,20,30])
        ax3[i,j].set_yticks([-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120])
        ax3[i,j].set_xlim([0,30])
        ax3[i,j].set_ylim([-120,120])
        ax3[i,j].set_facecolor('#DEE1E6')
        ax3[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax3[i,j].legend(loc=2)
        ax3[i,j].annotate('%s)' % subplotnum[count-1], xy=(26,97),fontsize=14)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig3.tight_layout()    
#fig3.savefig('MDOF_Relative_Displacement.png', bbox_inches= 'tight',dpi=300)    
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig3.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure14.tiff', bbox_inches= 'tight',dpi=300)

fig4, ax4 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        pbrotE = np.argmax(abs(baserot['%d' % count]))
        pbrotN = np.argmax(abs(sdofr['%d' % count]))
        
        ax4[i,j].plot(timearray['%d' % count], baserot['%d' % count], alpha= 1, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax4[i,j].plot(sdoft['%d' % count], sdofr['%d' % count], alpha= 0.75, color= '#FF8787', lw= 1.5, label= 'Numerical model')
        ax4[i,j].plot((0,30),(0.0027,0.0027),ls='--',color='k',lw=0.5,label=r'$0.01 \alpha$')
        ax4[i,j].plot((0,30),(-0.0027,-0.0027),ls='--',color='k',lw=0.5)
        #ax4[i,j].plot(timearray['%d' % count][pbrotE],baserot['%d' % count][pbrotE],'o',mec='limegreen',mfc='none')
        #ax4[i,j].plot(sdoft['%d' % count][pbrotN],sdofr['%d' % count][pbrotN],'o',mec='#FF8787',mfc='none')
        
        
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax4[i,j].set_xlabel('Time (s)')
        if count==1 or count==6 or count==11:
            ax4[i,j].set_ylabel('Base rotation (rad)')
        ax4[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, -0.06), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax4[i,j].annotate('Relative error of peak = '+str(round((abs(abs(sdofr['%d' % count][pbrotN])-abs(baserot['%d' % count][pbrotE]))/abs(baserot['%d' % count][pbrotE]))*100,1))+'%',
                      xy = (1, -0.07), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax4[i,j].set_xticks([0,10,20,30])
        ax4[i,j].set_yticks([-0.08,-0.06,-0.04,-0.02,0,0.02,0.04,0.06,0.08])
        ax4[i,j].set_xlim([0,30])
        ax4[i,j].set_ylim([-0.08,0.08])
        ax4[i,j].set_facecolor('#DEE1E6')
        ax4[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax4[i,j].legend(loc=2)
        ax4[i,j].annotate('%s)' % subplotnum[count-1], xy=(26,97),fontsize=14)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig4.tight_layout()    
#fig3.savefig('MDOF_Relative_Displacement.png', bbox_inches= 'tight',dpi=300)    
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig4.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/FigureBR.tiff', bbox_inches= 'tight',dpi=300)



'''

ratioS2 = [0.47278745856888066, 0.4155722460407499, 0.4699271860653164, 0.4596888186593633, 0.4494619542413561,
0.4451223127522226, 0.43744807274700975, 0.9068209956231483, 0.902647531869159, 0.8080042858262924, 0.9435616494250596,
0.9195416530693364, 1.1276079462263697, 0.8489606469997072, 0.8509329899403957]

ratioS2V = [1.0430196018065219, 0.9871337094328179, 1.0545155828236716, 1.049294113488027, 1.0680129774305556, 1.0715499752275157,
1.0549212705647557, 1.1734893322861075, 1.1018691346789102, 1.3393309107409137, 1.3554068613202979, 1.235237091043132,
1.4517115932988771, 1.068638851558815, 1.255403077028685,]

np.mean(ratioS2V)

testnum3 = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]

len(ratioS2)

dispmaxE = np.zeros(len(testnumber))
dispmaxN = np.zeros(len(testnumber))

VmaxE = np.zeros(len(testnumber))
VmaxN = np.zeros(len(testnumber))

for i in range(len(testnumber)):
    dispmaxE[i] = max(abs(dispdacm['%d' % (i+1)]))
    dispmaxN[i] = max(abs(sdofd['%d' % (i+1)]))
    VmaxE[i] = max(abs(accelcm['%d' % (i+1)]*m))
    VmaxN[i] = max(abs(sdoff['%d' % (i+1)]))
    print(abs(dispmaxN[i]-dispmaxE[i])/dispmaxE[i])
   



fig6, ax6 = plt.subplots(1,2, figsize= (20,8))
ax6[0].plot(testnum3,dispmaxE/dispmaxN,'o',color='#2D708EFF',label='Specimen 1')
ax6[0].plot(testnum3,ratioS2,'o',color='#440154FF',label='Specimen 2')
ax6[0].plot((0,16),(1,1),ls='--',color='k')
ax6[0].set_ylim(0.25,1.4)
ax6[0].set_xlim(0,16)
ax6[0].set_xticks(testnum3)
ax6[0].set_facecolor('#DEE2E6')
ax6[0].set_xticklabels(['2','3','4','5','6','7','8','9','11','12','13','14','15','17','18'],fontsize=14)
ax6[0].set_ylabel(r'$PCRD_{exp}/PCRD_{num}$',fontsize=14)
ax6[0].set_xlabel('Test',fontsize=14)
ax6[0].grid(color= 'black', ls= '--', lw= 0.5)
ax6[0].legend(loc=2,fontsize=14)
ax6[0].annotate('a)', xy=(15.2,0.3),fontsize=16)

ax6[1].plot(testnum3,VmaxE/VmaxN,'o',color='#2D708EFF',label='Specimen 1')
ax6[1].plot(testnum3,ratioS2V,'o',color='#440154FF',label='Specimen 2')
ax6[1].plot((0,16),(1,1),ls='--',color='k')
ax6[1].set_ylim(0.8,2.4)
ax6[1].set_xlim(0,16)
ax6[1].set_xticks(testnum3)
ax6[1].set_facecolor('#DEE2E6')
ax6[1].set_xticklabels(['2','3','4','5','6','7','8','9','11','12','13','14','15','17','18'],fontsize=14)
ax6[1].set_ylabel(r'$Vb_{exp}/Vb_{num}$',fontsize=14)
ax6[1].set_xlabel('Test',fontsize=14)
ax6[1].grid(color= 'black', ls= '--', lw= 0.5)
ax6[1].legend(loc=2,fontsize=14)
ax6[1].annotate('b)', xy=(15.2,0.85),fontsize=16)

fig6.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure18.tiff',dpi=300)
'''
Ta = np.linspace(0,1,100)

SaE = {}
SdE = {}
SaN = {}
SdN = {}

for i in range(len(testnumber)):
    SaE['%d' % (i+1)],SdE['%d' % (i+1)] = spectrumAD(timecmc['%d' % (i+1)],accelcm['%d' % (i+1)],Ta)
    SaN['%d' % (i+1)],SdN['%d' % (i+1)] = spectrumAD(sdoft['%d' % (i+1)],sdofa['%d' % (i+1)]/1000,Ta)


fig5, ax5 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr):
    for j in range(nc):
        ns= nc*i+(j+1)
        relES = np.mean(abs(SaN['%d' % ns][30:]-SaE['%d' % ns][30:])/SaE['%d' % ns][30:]*100)
        ax5[i,j].set_facecolor('#DEE1E6')
        ax5[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        #ax4[i,j].set_xticks([0,10,20,30])
        #ax4[i,j].set_yticks([0,100,200,300,400,500,600,700,800])
        #ax4[i,j].set_xlim([0,30])
        #ax4[i,j].set_ylim([0,800])
        ax5[i,j].plot(Ta,SaE['%d' % ns],color='limegreen',lw=1.5,label='Experiment')
        ax5[i,j].plot(Ta, SaN['%d' % ns],color='#FF8787', lw= 1.5, label= 'Numerical model')
        ax5[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[ns-1],intensity[ns-1]),
                     xy = (0.02, location4[ns-1]), 
                     xycoords = 'data',
                     xytext = (+0, +0), 
                     textcoords = 'offset points', 
                     fontsize = 10,
                     # weight= 'bold',
                     color = 'black',
                     )
        ax5[i,j].annotate('Mean relative error = '+str(round(relES,1))+'%',
                     xy = (0.02, location4[ns-1]-0.5), 
                     xycoords = 'data',
                     xytext = (+0, +0), 
                     textcoords = 'offset points', 
                     fontsize = 10,
                     # weight= 'bold',
                     color = 'black',
                     )
        ax5[i,j].annotate('%s)' % subplotnum[ns-1], xy=(0.9,5.5),fontsize=14)
        if ns == 1 or ns == 6 or ns==11:
            ax5[i,j].set_ylabel('Spectral Acceleration (g)')
        if ns == 11 or ns == 12 or ns==13 or ns==14 or ns==15:
            ax5[i,j].set_xlabel('Period (s)')
        if ns == 1:
            ax5[i,j].legend(loc=2)
        
        if ns==13:
            ax5[i,j].set_xticks([0,0.2,0.4,0.6,0.8,1])
            ax5[i,j].set_yticks([0,1,2,3,4,5,6])
            ax5[i,j].set_xlim([0,1])
            ax5[i,j].set_ylim([0,6])
        else:
            ax5[i,j].set_xticks([0,0.2,0.4,0.6,0.8,1])
            ax5[i,j].set_yticks([0,1,2,3,4,5,6])
            ax5[i,j].set_xlim([0,1])
            ax5[i,j].set_ylim([0,6])


fig5.tight_layout()   

#plt.savefig('MDOF_InSpectraAcc.tiff',dpi=300)
fig5.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure16.tiff', bbox_inches= 'tight',dpi=300)
